Solids and supersolids of three-body interacting polar molecules in an optical lattice 
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We study the physics of cold polar molecules loaded into an optical lattice in the regime of strong three-body 
interactions, as put forward recently by Buchler et al. [Nat. Phys. 3, 726 (2007)]. To this end quantum Monte 
Carlo simulations, exact diagonalization and a semiclassical approach are used to explore hardcore bosons 
on the 2d square lattice which interact solely by long ranged three-body terms. The resulting phase diagram 
shows a sequence of solid and supersolid phases. Our findings are directly relevant for future experimental 
implementations and open a new route towards the discovery of a lattice supersolid phase in experiment. 
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Introduction Strongly correlated systems studied in con- 
densed matter physics or in atomic physics are usually domi- 
nated by two-body interactions. The paradigm models are the 
fermionic and the bosonic Hubbard model which include a 
local two-body density-density interaction or the Heisenberg 
model consisting of two-body spin exchanges. These standard 
models are able to describe an enormous number of physical 
phenomena since the simultaneous interaction between more 
than two particles is small in most cases because it arises only 
in higher order of perturbation theory. Nevertheless multi- 
body interactions are present and can have profound effects on 
the physics of a system, e.g. ring-exchange processes being 
responsible for the rich nuclear magnetism of Helium 3 [1], 
the accurate description of undoped high-T c compounds and 
cuprate ladders requires four-spin interactions [2, 3, 4, 5], 
while three-body (3B) exchanges appear naturally in the con- 
text of two atomic species in a frustrated optical lattice topol- 
ogy [6]. 

On the theoretical side the study of microscopic models 
with multi-particle interactions is a very active and fruitful 
line of research. Exotic quantum ground states and decon- 
fined criticality [7] can possibly be triggered by such interac- 
tions [8, 9, 10]. Furthermore one can expect fractionalization 
of elementary excitations, e.g. spin liquid states in quantum 
magnets [11], topological ordered states as discussed in the 
context of quantum computation [12] or fractional quantum 
Hall states [13, 14, 15]. The major obstacle on the way to- 
wards an experimental confirmation of these fascinating pre- 
dictions is usually the requirement of dominating multi-body 
interactions, which is hard to achieve in a condensed matter 
setting. The field of ultracold gases loaded into optical lat- 
tices opens now a new perspective to overcome these difficul- 
ties. It has recently been shown that ultracold gases of polar 
molecules confined to optical lattices can be tuned to a regime 
where the interactions are solely of 3B type [16]. In contrast 
to conventional Hubbard or Heisenberg models having mostly 
short-range interactions, the 3B density interactions put for- 
ward in [16] decay only slowly in space, owing to the un- 
derlying dipolar interactions. It is an important task to study 
systems based on these novel multi-body interactions and to 




FIG. 1: (Color online) Illustration of the leading 3B interactions 
Wijk with amplitude larger than 0.2, defining the minimal model 
studied here. 

uncover the nature of unconventional phases harbored in their 
phase diagram. 

In the present Letter, we achieve a step in this direction by ex- 
ploring the potentially most relevant case of hardcore bosons 
on a square lattice interacting only through 3B forces. We per- 
form a comprehensive numerical study, based on a semiclas- 
sical approximation (SCA), exact diagonalization (ED) and 
quantum Monte Carlo simulations (QMC) to derive the re- 
sulting zero temperature phase diagram. We reveal a rich se- 
quence of solids, supersolids and phase separation as the den- 
sity n is tuned from to 1. Interestingly we find a stable, 
extended checkerboard supersolid (CSS) phase around den- 
sity n = 1/2. Such lattice supersolids are currently a topic 
of great interest both in the fields of cold atoms and quantum 
magnetism (see e.g. Refs. 17, 18). 

Model We consider hard-core bosons hopping on the two- 
dimensional square lattice including 3B density interactions 
as put forward in [16]: 

H = -t^T(b\bj + h.c.) - fi^Tni + - ^ U,,;. ''•''.;"< 

(1) 

where rii — b i b i is the boson density at site i, /i is the chemi- 
cal potential, t is the nearest-neighbor hopping amplitude, and 
W^k labels the 3B interactions. 

The 3B interactions derive from the dipolar forces between the 
polar molecules under the additional influence of microwave 



2 




FIG. 2: (Color online) Central plot: Phase diagram of the minimal model obtained within the semiclassical approach as a function of the 
chemical potential n and the nearest-neighbor hopping amplitude t. Grey (white) regions are superfluid (phase separated), and light (dark) 
blue denotes supersolids (solids). Surrounding plots: schematic representation of the nature of some of the solid and supersolid phases. The 
greyscale of the circles represents the filling (white = empty, black = full), while the length and the direction of the red arrows denotes the 
amplitude and the phase of the superfluid component. The blue lines highlight the unit cell of the different structures. 



fields [16], and retain some of its character, especially the long 
range nature. The general expression for the amplitudes Wijk 
is given by: 



W m = W 



1 

IR, - R^lRt - R k \ 



permutations 



(2) 

The coefficient Wo depends on the microscopic setup and is 
discussed in Ref. 16. In the following the energy scale is set 
by the 3B interactions and we thus put Wo = 1. Note that 
the spatial dependence of the interactions is such that the re- 
pulsion between 3 particles with a mutual distance of order 
R amounts only to 1/i? 6 . If however two particles are close, 
while the third is at distance R, then the interactions only de- 
cay as 1/i? 3 , resembling the decay of the underlying dipo- 
lar interactions. Based on these considerations we expect the 
3B interactions to have a stronger effect at high densities than 
at very low density. Furthermore by applying a particle-hole 
transformation in the regime of densities close to n = 1 one 
can effectively map the problem to a low-density, two-body 
dipolar gas of holes. So the most challenging regime remains 
for densities n ~ 1/2, where the full structure of the 3B inter- 
actions is important. 

The plan of the paper is to study first a minimal model where 
the range of the 3B interactions is limited to the 5 types of 
terms illustrated in Fig. 1 (i.e. Wijk > 0-2). We map out 
the phase diagram using a semiclassical approach, and con- 
firm the main findings for selected parameters by numerical 
ED and QMC simulations. Then we corroborate the utility 
of the minimal model by including all terms with amplitudes 
W^k > 10~ 3 in the semiclassical approach. 

Semiclassical approximation ( SCA ) The SCA maps Eq. 1 
to a spin 1/2 model using the exact Matsubara-Matsuda [19] 
representation of hardcore bosons S + = b, S~ = b\ and 
S z = 1/2 — n. The resulting spin Hamiltonian is studied in 
the classical limit by replacing the quantum spins by classical 



vectors of length 1/2 on a sphere. The classical ground state 
is obtained in the thermodynamic limit by numerically deter- 
mining the global energy minimum among all non-equivalent 
unit cells with up to 32 sites. The spin structures minimizing 
the energy are mapped back to the boson problem and cor- 
respond typically to superfluid, solid or supersolid phases of 
varying spatial complexity. A supersolid is a phase breaking 
simultaneously the U(l) gauge symmetry (superfluid) and the 
underlying translational symmetry of the lattice (solid). The 
method is computationally much less expensive than the ED 
or QMC simulations, and can therefore be used to efficiently 
map out the phase diagram. 

The resulting SCA phase diagram for the minimal model as a 
function of t and p, is shown in Fig. 2. For large t, the system 
corresponds to basically non-interacting hardcore bosons and 
is thus expected to be superfluid for all densities. In the op- 
posite limit t = only commensurate solid phases are found. 
The density n(p) displays a simple series of plateaux which 
are separated by first order transitions, i.e. jumps in the den- 
sity. We find plateaux at n = 1/3, 1/2, 5/8, 2/3, 3/4, 4/5, 
5/6, and 7/8 for the minimal model. The much richer struc- 
ture above n = 1/2 compared to low densities is a conse- 
quence of the particle-hole asymmetry discussed above. The 
specific structure of some of the plateaux are illustrated in 
Fig. 2. Due to the finite range of the truncated interactions 
the plateaux at 5/8,5/6 and 7/8 exhibit a residual degener- 
acy in the limit t = 0, which is expected to be lifted either by 
the longer range couplings (see below) or an order by disor- 
der mechanism driven by the quantum fluctuations at finite t, 
which is however beyond the present SCA approach. Within 
the SCA the quantum melting of the various solids takes place 
for values t < 2. The physics below n = 1/2 consists of a su- 
perfluid developing from low densities as t increases, a solid 
at n = 1/3, which is destroyed by a rather small amount of 
hopping, as well as a puzzling supersolid without correspond- 
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ing n = 3/8 plateau [18]. Finally the checkerboard solid (CS) 
at n = 1/2 extends to densities below 1/2 by forming a stable 
CSS with a non-trivial dependence on t, leading to a maximal 
region of stability around t w 0.6. For values t > 0.8 the tran- 
sition from the superfiuid below to the CS is direct and first or- 
der. The physics above the CS is even richer. Above n = 1/2, 
we find a CSS in a large range of t values. The plateau at 5/8 
also has a corresponding supersolid for a density range below 
5/8, but this supersolid is very compressible. Furthermore we 
find supersolids for densities just below n = 2/3 and 3/4 in a 
small window of t. 

The SCA phase diagram of the minimal model is very rich, 
including several supersolids of different spatial struture. We 
now proceed to a numerically exact treatment of the minimal 
model to confirm the main findings, i.e. the basic solids and 
the ft = 1/2 supersolids. 

Numerical simulations In the following we use QMC and 
ED in order to corroborate the predictions made by the SCA. 
We focus first on the case t = 0.5 displaying most of the fea- 
tures in the SCA calculation and briefly comment on results 
obtained for other t values. The ED calculations were per- 
formed on square clusters up to 36 sites. The QMC simula- 
tions are based on a modified [21] stochastic series expansion 
(SSE) [20, 22] code and the ALPS libraries [23]. We restrict 
ourselves to the density range < n < 5/8 for t = 0.5. 
On the one hand because the SSE algorithm based on the di- 
rected loop update is not particularly efficient in exploring the 
solid phases with large unit cells found at higher densities, and 
on the other hand it is also the most interesting density range 
since it displays a sizable CSS phase in the SCA approach. 
The lowest temperature was typically T = 0.05, which is rep- 
resentative of the ground state for the shown quantities. Sys- 
tems sizes went up to N = 12 x 12 = 144 sites. 

The numerical data for the density n as function of p, is 
shown in the top panel of Fig. 3. The systems starts to fill 
at /ii = —At = —2 and the density behaves smoothly up to 
p w 4.2 where the system phase separates between a super- 
fluid component at density n « 0.4 and the CS at n = 1/2. 
Note that there is no n = 1/3 plateau present at t = 0.5. 
The n = 1/3 plateau is recovered for t = 0.25 (see below) 
implying that at least some differences to the SCA can be re- 
solved in terms of a renormalized t. At // « 13.6 the CS gets 
doped with particles which condense without destroying the 
solid and thus form a CSS. The supersolid remains stable up 
to/is; 16.3(1) where phase separation occurs anew, this time 
between the CSS and the n = 5/8 solid. 

The phases just discussed are determined by measurements 
of the superfiuid stiffness p s — 2 ^ L2 (W% 



Wy) where 



W x and W y are the total winding numbers in x and y direc- 
tions, the condensate fraction po = lim^oo (&J& ) as well 
as the checkerboard charge order parameter S(tt,tt)/N = 
-^2 2wi j 'X~ ( n i n j) displayed in the central and bottom 
panels of Fig. 3. At low densities < n < 0.4 the system 
is indeed superfiuid, with a finite superfiuid stiffness p s and 
condensate fraction pq. At very low densities both quanti- 
ties are in very good quantitative agreement with the values 
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FIG. 3: Upper graph: density n as a function of the chemical poten- 
tial \i for the minimal model (c.f. Fig. 2) at t = 0.5 up to density 
n — 5/8. Lower graph, upper panel: checkerboard order parameter 
S(tv, tt)/N from ED (N = 32) and QMC (N = 12 X 12) simula- 
tions. Lower panel: superfiuid stiffness p s and condensate fraction 
po obtained by QMC. Inset: p s at fixed j3 — 20 and p = 15 for 
L = 6,8, ...,16. 



obtained for noninteracting hardcore bosons [24]. At density 
n = 1/2 the checkerboard order is highlighted by the struc- 
ture factor data obtained with ED and QMC for different sys- 
tem sizes. As one dopes the CS with additional particles a 
sizable CSS emerges for 1/2 < n < 0.59, therefore confirm- 
ing the SCA prediction. Finite size effects are small in the 
CSS phase, as witnessed by S(tt, tt) /N being essentially un- 
changed from N = 32 (ED) to N = 144 (QMC), while the 
finite size extrapolation of the stiffness p s (QMC, inset) con- 
verges to a finite value, showing that the CSS is stable in the 
thermodynamic limit. 

The stability of the CSS is surprising, since it is commonly 
believed that the CSS is unstable for hardcore bosons with 
only nearest-neigbor hopping [25]. In the present case it can 
however be shown that the instability towards domain-wall 
formation [26] is absent, therefore providing an explanation 
for the stability of the CSS. We remark that in the CSS phases 
the solid order is more pronounced than the superfiuid compo- 
nent, due to the vicinity of the solid. We therefore expect the 
supersolid to first give way to a solid phase by a Kosterlitz- 
Thouless transition, followed by an Ising transition to a nor- 
mal bose liquid with increasing temperature [18]. 

We have performed ED simulations for t = 0.25, 0.75 and 
1 to further check for the presence of phases predicted by the 
SCA [27]. At t = 0.25 we find evidence for a n = 1/3 
plateaux, a CSS below as well as above the n = 1/2 CS. At 
t = 0.75 the CSS above the n = 1/2 plateau and the phase 
separation below the solid are reduced in density extent, com- 
pared to t — 0.5. Finally at t = 1, using QMC and ED, we 
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FIG. 4: (Color online) Phase diagram obtained by SCA for different 
values of t including all 3B terms with amplitudes Wijk > 0.001. 
Grey regions are superfluid, and light (dark) blue denotes supersolids 
(solids), which were already present in the minimal model c.f. Fig. 2. 
The white regions encompass phase separation and further solids and 
supersolids. 



find an— 1/2 and n — 2/3 plateau, while the CSS above 
n = 1/2 is tiny if present at all. The remaining regions are 
superfluid. The phase transitions between the superfluid and 
the solids at t — 1 deserve further study. 
Although the SCA approach does not treat the quantum fluc- 
tuations quantitatively correctly, our numerical investigations 
confirms that many of the qualitative SCA phases are correct 
and indeed present. Based on this validation we now address 
the effect of the finite range approximation. 

Effect of longer range interactions The numerical solu- 
tion of the Hamiltonian (1) including the 3B interactions at all 
distances is a formidable task. Even the classical problem for 
t = is non-trivial due to the discrete structure imposed by 
the square lattice. Here we are merely interested whether the 
generic features found in the truncated model are stable upon 
the inclusion of the long-range nature of the 3B interactions. 
To this end we use the SCA including all 3B interactions in the 
Hamiltonian with amplitudes W^k > 0.001. Minimizations 
are done on all clusters up to 24 sites. While these clusters 
might still be too small to represent all structures at small t, 
they are amply sufficient to confirm that many of the solids 
and supersolids present in the minimal model survive the in- 
clusion of the long-range 3B couplings. The phase diagram is 
shown in Fig. 4. Most importantly we confirm the stability of 
all plateaux of the minimal model, plus the supersolids below 
and above the CS, as well as supersolids below the n = 2/3 
and 3/4 solids. While in some solids at large densities (e.g. 
n = 5/6 and 7/8) the degeneracy on the classical level is 
lifted, we find that other solids (n = 5/8 and 3/4) change 
slightly the charge order pattern. This reflects the general ten- 
dency of the 3B interactions at large densities to form a trian- 
gular Wigner crystal of holes. The optical lattice in the square 
geometry then acts as an incommensurate substrate, possibly 
giving rise to physics similar to the Frenkel-Kontorova model 
and to slow equilibration (glassiness) due to many almost de- 
generate charge configurations. 



Conclusion We studied a model of hardcore bosons on the 
square lattice interacting solely by slowly decaying 3B inter- 
actions, which is directly relevant for future experiments on 
ultracold gases of polar molecules confined to optical lattices. 
We find a rich phase diagram consisting of many solid, super- 
fluid, and supersolid phases. The long-range nature of the 3B 
interactions results in a zoo of crystalline phases in the limit of 
small t. The large number of competing states will probably 
also lead to difficulties in the equilibriation of the experimen- 
tal system in this limit. 

The most important finding of our work is that a sys- 
tem which only contains 3B interactions realizes supersolid 
phases. Extended supersolid phases exist around the n = 1/2 
CS. Our findings therefore suggest that ultracold gases of po- 
lar molecules confined to optical lattices are promising candi- 
dates to observe for the first time a supersolid on a lattice in 
experiment. 
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